Использовалась библиотека stats

Библиотеки

library(Rssa)
library(ggplot2)
library(mFilter)
data(USUnemployment)

Периодограммы

к 14.03

Эффект растекания синуса

Растекание — эффект, возникаемый, если длинна ряда не будет кратна периоду. Найдем длинну нашего ряда:

w <- 0.1  # Тут задаем частоту, которая соответствует периоду T = 10
x <- 1:100
y <- sin(2 * pi * w * x)
length(y)

## [1] 100

Получаем длинну ряда 100, построим периодограмму для ряда у которого длинна кратна периоду:

sp <- spec.pgram(y, spans = NULL, detrend = FALSE, log = "no", fast = FALSE, taper = 0)

Видим, что эффекта растекания на периодограмме нет, пик единственный и соответсвует заданной частоте. Теперь попробуем изменить длинну исходного ряда. Можно добавить несколько нулевых наблюдений в ряд, но лучше убрать уже имеющиеся, чтобы не искажать картину нулями:

yCut <- y[1:(length(y) - 3)] # Убираем три наблюдения
length(yCut)

## [1] 97

Теперь длинна ряда 97, что не кратно периоду T = 10, теперь посмотрим на его периодограмму:

sp <- spec.pgram(yCut, log = "no", detrend = FALSE, fast = FALSE, taper = 0)

Видим, что пик начал расползаться по соседним частотам.

Генерация белого шума

Шум — это случайная составляющая временного ряда. Если средний вклад всех частот одинаковый, то шум называется белым. Свойство периодограммы белого шума: значение периодограммы имеет экспоненциальное распределение с одним и тем же средним.

n <- 1000
wnoise <- rnorm(n, 0, 1)
spec.pgram(wnoise, detrend = FALSE, log = "no", fast = FALSE, pad = FALSE,taper = 0)

Так как для шума вообще нет детерминирующих частот, мы видим из периодограммы, что каждая частота входит примерно с одинаковой силой. С увеличением количества наблюдений будет увеличиваться значение каждой частоты на сетке:

nLarge <- 10000
wnoiseLarge <- rnorm(nLarge, 0, 1)
spec.pgram(wnoiseLarge, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE, taper = 0)

Генерация красного шума

Отличие красного шума от белого в том, что предыдущие значения влияют на следующие, т.е. имеют ненулевую корреляцию, в отличие от белого шума.

\(r_n = cor \cdot r_{n-1}+\sqrt{(1-cor^2)} \cdot \omega\)

w0 <- wnoise[1]
wnoise <- wnoise[2:n]
cor <- 0.8
rnoise = Reduce(function(prev_v, next_v) cor * prev_v + next_v * sqrt(1 - cor^2), wnoise, w0, accumulate = T)

spec.pgram(rnoise, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

Построение периодограммы для реальных данных

График безработицы в США

О данных: https://link.springer.com/chapter/10.1007/978-1-4612-5098-2_66 408 ежемесячных наблюдений безработицы в США среди мужчин Длину окна берем равную 12 в соотвествии с годовой периодичностью. Посмотрим на их график:

library(Rssa)
data(USUnemployment)

unempl.male <- USUnemployment[, "MALE"]

plot(unempl.male, type = "l")

Периодограмма безработицы в США

По графику можно сказать, что временной ряд достаточно волатильный, имеет сложный тренд, по периодограмме можно сделать выводы о частотах входящих в него:

spec.pgram(unempl.male, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

Видим, что наибольший вклад имеют низкие частоты, которые мы относим к шуму, также есть пики в точках w = 1/12 и w = 2/12, что говорит нам о наличии годовой и полугодовой сезонности.

Фильтры

к 14.03

Функция для фильтра скользящего среднего

movingAverage <- function(x, n=1, centered=TRUE) {
      
      if (centered) {
            before <- floor  ((n-1)/2)
            after  <- ceiling((n-1)/2)
      } else {
            before <- n-1
            after  <- 0
      }

      s     <- rep(0, length(x))
      count <- rep(0, length(x))
      
      new <- x
      count <- count + !is.na(new)
      new[is.na(new)] <- 0
      s <- s + new

      i <- 1
      while (i <= before) {
            new   <- c(rep(NA, i), x[1:(length(x)-i)])
            
            count <- count + !is.na(new)
            new[is.na(new)] <- 0
            s <- s + new
            
            i <- i+1
      }
      
      i <- 1
      while (i <= after) {
            new   <- c(x[(i+1):length(x)], rep(NA, i))
            
            count <- count + !is.na(new)
            new[is.na(new)] <- 0
            s <- s + new
            
            i <- i+1
      }
      s/count
}

АЧХ

afc <- function(filter, omega)
{
      k <- seq_along(filter) - 1
      h <- function(o) sum(rev(filter)*exp(-k*1i*o))
      abs(sapply(omega, h))
}

Применение фильтра скользящего среднего к модельным данным

Задаем модель

n = 100
x <- 1:n
period <- 10

modelTs <- (2*x + 3) +
 5*cos(2*pi*(1/period)*x) + 5*cos(2*pi*(2/period)*x) + 5*cos(2*pi*(3/period)*x) + 5*cos(2*pi*(4/period)*x) +
  + rnorm(n)

plot(modelTs, type = "l",  main = "Исходный временной ряд")

Исходная периодограмма

sp <- spec.pgram(modelTs, detrend = FALSE, log = "no", fast = FALSE, pad = FALSE, taper = 0)

Применяем фильтр

plot(modelTs, main = "Исходный ряд + Cкользящее среднее", type = "l")
legend(x = "topleft",                                # Position
       legend = c("timeseries", "MA(5)", "MA(10)"),  # Legend texts
       lty = c(1, 1, 1),                             # Line types
       col = c("black", "red", "blue"),              # Line colors
       lwd = 2)                                      # Line width

maHalf <- movingAverage(modelTs, 5, TRUE)
lines(maHalf, col = "red")

maFull <- movingAverage(modelTs, 10, TRUE)
lines(maFull, col = "blue")

Периодограмма MA(5)

sp <- spec.pgram(maHalf ,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

sp <- spec.pgram(maHalf ,detrend = FALSE, log = "yes",fast = FALSE, pad = FALSE,taper = 0)

Периодограмма MA(10)

sp <- spec.pgram(maFull ,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

sp <- spec.pgram(maFull ,detrend = FALSE, log = "yes",fast = FALSE, pad = FALSE,taper = 0)

АЧХ фильтра скользящего среднего

freq <- seq(0, pi, 0.001)
filt <- rep(1/period, period)
omega <- freq/2/pi

str <- c("0", "1/10", "2/10", "3/10", "4/10","5/10")
plot(afc(filt, freq) ~ omega, type = "l", xlab = "Frequency", ylab = "Frequency response", xaxt = "n")
title(main = "АЧХ фильтра скользящего среднего MA(10)")
axis(1, at = seq(0, 0.5, by = 1/10), las = 2, labels = sprintf("%s", str))

period = 5
filt <- rep(1/period, period)

plot(afc(filt, freq) ~ omega, type = "l", xlab = "Frequency", ylab = "Frequency response")
title(main = "АЧХ фильтра скользящего среднего MA(5)")

Применение фильтра к реальным данным

График безработицы в США Периодограмма безработицы в США

Для выделения тренда выставим длинну окна равную периоду ряда (T = 12). Наложим график скользящего среднего на исходный:

plot(unempl.male, main = "Исходный ряд + Cкользящее среднее", type = "l")
legend(x = "topleft",                                
       legend = c("timeseries", "MA(12)", "MA(24)"),
       lty = c(1, 1, 1, 1),                           
       col = c("black", "blue", "red"),    
       lwd = 2)                                     

ma12 <- movingAverage(unempl.male, 12, TRUE)
ma24 <- movingAverage(unempl.male, 24, TRUE)
lines(ma12, col = "blue")
lines(ma24, col = "red")

Теперь построим периодограмму скользящего среднего:

spec.pgram(ma12, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)

SSA

к 14.03

Выделение тренда с помощью SSA

Применим SSA к реальным данным, ширину окна можно взять L = K/2 и кратную периоду (условие ассимптотической разделимости). Построим матрицу взвешенных корреляций, чтобы проанализировать компоненты:

Матрица взвешенных корреляций

tsSSA <- ssa(unempl.male, L = length(unempl.male) %/% 2)
tsSSACor <- wcor(tsSSA, groups = 1:30)
plot(tsSSACor)

Видим, что после 4 компоненты начинается зашумленность, при этом 15-я компонента достаточно хорошо отделяется от соседних.

Собственные вектора

plot(tsSSA, type="vectors", idx = 1:20)

По парам собственных векторов можно сказать о векторах периодичностей:

plot(tsSSA, type="paired", idx = 2:20)

Видим правильные фигуры в компонентах (5, 6) и (12, 13):

tsSSASeason <- reconstruct(tsSSA, groups = list(season = c(5, 6, 12, 13)))
spec.pgram(tsSSASeason$season, log = 'no', fast = FALSE, taper = 0, detrend = FALSE)

Видим, что выделились пики соответсвующие годовой и полугодовой периодичности, как на исходной периодограмме.

tsSSASignal <- reconstruct(tsSSA, groups = list(c(1:4, 7:10), c(5, 6, 12, 13)))

plot(tsSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

FOSSA

Алгоритм улучшает сильную разделимость.

FOSSA: матрица корреляций

tsFOSSA <- fossa(tsSSA, nested.groups = c(1:13))
tsFOSSACor <- wcor(tsFOSSA, groups = 1:30)

plot(tsFOSSACor)

FOSSA: собственные вектора

plot(tsFOSSA, type = "vectors", idx=1:15)

FOSSA: итоговое разложение

tsFOSSASignal = reconstruct(tsFOSSA, groups <- list(c(5:13), 1:4))

plot(tsFOSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

Корреляции SSA vs. FOSSA

plot(tsSSACor)

plot(tsFOSSACor)

Сезонности SSA vs. FOSSA

matplot(data.frame(tsSSASignal$F2, tsFOSSASignal$F2), type = 'l', col=c("red","green"), lty=c(1,1))

Тренд SSA vs. FOSSA

matplot(data.frame(c(unempl.male), tsSSASignal$F1, tsFOSSASignal$F1), type = 'l', col=c("black", "red","green"), lty=c(1,1))

Различные способы выделения тренда

к 21.03, 28.03

Тренд по POLY

trend.poly3 = lm(unempl.male ~ poly(1:408, degree = 3))
trend.poly5 = lm(unempl.male ~ poly(1:408, degree = 5))
trend.poly7 = lm(unempl.male ~ poly(1:408, degree = 7))


plot(c(unempl.male), main = "Исходный ряд + POLY", type = "l")
legend(x = "topleft",                                
       legend = c("timeseries", "POLY(3)", "POLY(5)", "POLY(7)"),
       lty = c(1, 1, 1, 1),                           
       col = c("black", "red", "blue", "green"),    
       lwd = 2)    

lines(trend.poly3$fitted.values, col="red", type = "l")
lines(trend.poly5$fitted.values, col="blue", type = "l")
lines(trend.poly7$fitted.values, col="green", type = "l")

Тренд по LOESS

LOESS метод основан на сглаживании с помощью построения локальных взвешенных линейных регрессий

index <- 1:408

fit.loess <- loess(unempl.male ~ index, span = 0.2, degree = 1)      
trend.loess <- predict(fit.loess)

plot(c(unempl.male), main = "Исходный ряд + LOESS", type = "l")
legend(x = "topleft",                                
       legend = c("timeseries", "LOESS"),
       lty = c(1, 1, 1, 1),                           
       col = c("black", "blue"),    
       lwd = 2)                                     
lines(trend.loess, col = "blue", type = "l")

Тренд по HP

chest = hpfilter(unempl.male, freq=60000, type="lambda")
trend.hp = chest$trend
plot(unempl.male, type="l")
lines(trend.hp, col="red")
legend("topleft", legend = c("Real TS","HP"), col = c("black","red"), lty=1)

Итоговый график со всеми способами выделения тренда

plot(unempl.male, type="l")
lines(trend.hp, col="red")
lines(ma24, col="blue")
lines(ts(trend.loess, start=c(1948, 1), frequency=12), col="orange")
lines(ts(trend.poly7$fitted.values, start=c(1948, 1), frequency=12), col="green")
legend("topleft", legend = c("Real TS","HP", "MA(24)", "LOESS", "POLY(5)"), col = c("black","red", "blue", "orange", "green"), lty=1)

Сравнение методов разложения

04.04

STL

STL имеет следующие параметры:

(inner) – число итераций внутреннего цикла (р).
(outer) – число итераций внешнего цикла. У нас в данных нет аутлаеров, поэтому outer = 0
(l.window) – сглаживающий параметр для low-pass фильтра. Возьмем равным 13. Рекомендуют брать ближайший нечетный к периоду.
(t.window) – сглаживающий параметр для тренда. Рекомендуют брать ближайшим нечетным к (1.5*period) / (1-(1.5/s.window)). Равен 23.

ns - (s.window) – сглаживающий параметр сезонности. Берем 13, так как годовая периодичность и параметр должен быть нечетным

tsSTL <- stl(unempl.male, s.window = 13, l.window = 13, outer = 0, inner = 1, t.window = 23)
plot(tsSTL)

Classical seasonal decomposition

tsCSD <- ssa(unempl.male, L = length(unempl.male) %/% 2, force.decompose =  FALSE, svd.method = "nutrlan")
fit.dec <- decompose(tsCSD, neig = 40)
tsCSDCor <- wcor(fit.dec, groups= 1:40)
plot(tsCSDCor)

plot(fit.dec, type = "vectors", idx = 1:20)

tsCSDSignal <- reconstruct(fit.dec, groups = list(trend = c(1:4, 7:10, 15), season = c(5, 6, 12, 13)))
plot(tsCSDSignal)

Итоговый график трендов разных разложений

plot(unempl.male, type="l")
lines(tsSSASignal$F1, col="red")
lines(tsFOSSASignal$F1, col="blue")
lines(tsCSDSignal$trend, col="orange")
lines(ts(data.frame(tsSTL$time.series)$trend, start=c(1948, 1), frequency = 12), col="green")
legend("topleft", legend = c("Real TS","SSA", "FOSSA", "CSD", "STL"), col = c("black","red", "blue", "orange", "green"), lty=1)

Примеры применения Toeplitz SSA и SSA with projections

SSA with projections

Projection SSA применяется для улучшения разделимости линейного тренда и сезонности. Используем смоделированные данные с линейным трендом (Применение фильтра скользящего среднего к модельным данным):

x <- 1:100
modelTs <- (2*x + 3) +
 + 5*cos(2*pi*(1/25)*x) + cos(2*pi*(1/5)*x) + rnorm(50)

plot(modelTs, type="l")

# Применим SSA с двойным центрированием для лучшего отделения линейного тренда
modelTsSSA <- ssa(modelTs, L = 25, column.projector='centering', row.projector='centering')
plot(modelTsSSA, type = "vectors")

plot(modelTsSSA, type = "paired")

modelTsSSASignal <- reconstruct(modelTsSSA, groups = list(1:2))

plot(modelTs, type = "l")
lines(modelTsSSASignal$F1, col = "red")

Toeplitz SSA

Теплицев SSA применяется как улучшение basic-SSA для стационарных рядов.

x <- 1:100
modelTs <- 5*cos(2*pi*(1/25)*x) + 2*cos(2*pi*(1/5)*x) + rnorm(50)

plot(modelTs, type="l")

modelTsSSA <- ssa(modelTs, kind = "toeplitz-ssa")

plot(modelTsSSA, type = "vectors")

plot(modelTsSSA, type = "paired")

modelTsSSASignal <- reconstruct(modelTsSSA, groups = list(1:4))

plot(modelTsSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

plot(modelTs, type = "l")
lines(modelTsSSASignal$F1, col = "red")

Прогнозы: рекурентный и векторный

02.05

Рекурентный прогноз - rforecast

forcastR <- rforecast(tsSSA, groups = list(trend = c(1:4, 7:10), season = c(5, 6, 12, 13)), len = 12, only.new = F)

plot(unempl.male, type="l")
lines(data.frame(forcastR)$trend, col="red")
legend("topleft", legend = c("Real TS","RForcast trend"), col = c("black","red"), lty=1)

C отсечением периода

Применяем SSA к ряду с отсеченным периодом

tsCut = ts(unempl.male[1:(length(unempl.male)-12)], start=c(1948, 1), frequency=12)

tsCutSSA <- ssa(tsCut,  L = length(tsCut) %/% 2)
plot(wcor(tsCutSSA, groups = 1:30))

plot(wcor(tsCutSSA, groups = 1:15))

plot(tsCutSSA, type="vectors", idx = 1:20)

tsCutSSASignal = reconstruct(tsCutSSA, groups <- list(c(1:3, 6:10, 14:15), c(4,5,12,13)))

plot(tsCutSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))

forcastRCut <- rforecast(tsCutSSA, groups = list(trend = c(1:3, 6:10, 14:15), season = c(4,5,12,13)), len = 12, only.new = F)
plot(c(tsCut), type="l")
lines(c(data.frame(forcastRCut)$trend), col="red")
lines(c(data.frame(forcastR)$trend), col="blue")
legend("topleft", legend = c("Cut TS","RForcast cut trend", "RForcast full trend"), col = c("black","red", "blue"), lty=1)

Векторный прогноз - vforecast

forcastV <- vforecast(tsSSA, groups = list(trend = c(1:4, 7:10), season = c(5, 6, 12, 13)), len = 12, only.new = F)

plot(cbind(unempl.male, forcastV$trend), plot.type = "single", 
     col = c("black", "red"), ylab = NULL)

legend("topleft", legend = c("Real TS","VForcast full trend"), col = c("black","red"), lty=1)

C отсечением периода

forcastVCut <- vforecast(tsCutSSA, groups = list(trend = c(1:3, 6:10, 14:15), season = c(4,5,12,13)), len = 12, only.new = F)


plot(c(tsCut), type="l")
lines(c(data.frame(forcastVCut)$trend), col="red")
lines(c(data.frame(forcastV)$trend), col="blue")
legend("topleft", legend = c("Real TS","VForcast cut trend", "VForcast full trend"), col = c("black","red", "blue"), lty=1)

Восстановление ряда по ЛРФ

25.04/02.05

Ряд 1: Экспонента

library(Matrix)

Находим ранг матрицы для экспоненциального ряда:

N = 10
LRF = 3.1^(1:N)
rankMatrix(hankel(LRF))[1]

## [1] 1

Ожидаемо, он равен одному. Далее находим сигнальные корни ряда:

LRF.ssa = ssa(LRF, L=2, method = "svd")
l = lrr(LRF.ssa, groups = list(1))
plot(l)

roots_my = roots(l)
K = length(roots_my)
S_n = LRF[1:K]
lin_sys = data.frame(S_n)
vars = matrix(nrow = K, ncol = K)
for(i in 1:K){
  vars[i]=roots_my[i]^(1:K)
}
lm0 = lm(formula = S_n ~ 0 + ., data = data.frame(vars))
print(lm0)

## 
## Call:
## lm(formula = S_n ~ 0 + ., data = data.frame(vars))
## 
## Coefficients:
## vars  
##    1

Находим значения по линейной рекурентной формуле:

pred_S_n = lm0$coefficients[[1]]*(roots_my[1]^(1:N))

print(LRF)

##  [1]     3.1000     9.6100    29.7910    92.3521   286.2915   887.5037
##  [7]  2751.2614  8528.9104 26439.6222 81962.8287

S_n = LRF[1:N]
lin_sys = data.frame(S_n)
vars = matrix(nrow = N, ncol = K)
for(i in 1:K){
  vars[i]=roots_my[i]^(1:N)
}

lm0 = lm(formula = S_n ~ 0 + ., data = data.frame(vars))
print(lm0)

## 
## Call:
## lm(formula = S_n ~ 0 + ., data = data.frame(vars))
## 
## Coefficients:
## vars  
##    1

Сравниваем исходный ряд с полученным

pred_S_n = lm0$coefficients[[1]]*(roots_my[1]^(1:N))

print(LRF)

##  [1]     3.1000     9.6100    29.7910    92.3521   286.2915   887.5037
##  [7]  2751.2614  8528.9104 26439.6222 81962.8287

print(pred_S_n)

##  [1]     3.1000     9.6100    29.7910    92.3521   286.2915   887.5037
##  [7]  2751.2614  8528.9104 26439.6222 81962.8287

Ряд 2: Линейная функция

N = 10

LRF = 2*(1:N)+5

Находим ранг матрицы:

rankMatrix(hankel(LRF))[1]

## [1] 2

Линейная функция задается двумя компонентами.

LRF.ssa = ssa(LRF, L=3, method = "svd")
l <- lrr(LRF.ssa, groups = list(1:2))
roots_my <- roots(l)
print(l)

## [1] -1  2
## attr(,"class")
## [1] "lrr"

print(roots_my)

## [1] 1.0000002 0.9999998

Получаем два комплексно-сопряженных корня

print(2 * pi / Arg(roots_my))

## [1] Inf Inf

print(Mod(roots_my))

## [1] 1.0000002 0.9999998

roots_my2 <- rep(Re(mean(roots_my)), length(roots_my))
print(roots_my2)

## [1] 1 1

Находим коэффиценты исходного ряда

m = length(roots_my2)
s_n = LRF[1:N]
vars <- matrix(nrow = N, ncol = m)
for (i in 1:m) {
  vars[ ,i] <- (1:N) ^ (i-1) * roots_my2 ^ (1:N)
}

lm0 <- lm(s_n ~ 0 + ., data = data.frame(vars))
print(lm0)

## 
## Call:
## lm(formula = s_n ~ 0 + ., data = data.frame(vars))
## 
## Coefficients:
## X1  X2  
##  5   2

Сравниваем исходный ряд с полученным

pred_S_n = lm0$coefficients[[1]][1] * (vars[,1]) + lm0$coefficients[[2]][1] * (vars[,2])
print(LRF)

##  [1]  7  9 11 13 15 17 19 21 23 25

print(pred_S_n)

##  [1]  7  9 11 13 15 17 19 21 23 25

Восстановление ЛРФ реального временного ряда

Проанализируем собственные числа ряда

tsSSA$sigma

##  [1] 367075.857  74505.552  63304.591  35713.403  34728.609  34583.735
##  [7]  27721.976  22861.296  22230.219  18781.365  18634.256  18365.270
## [13]  18171.464  14245.491  13430.646  13113.963  10353.466  10239.538
## [19]   9456.080   9280.711   9009.749   8960.235   8234.105   8146.040
## [25]   7272.790   6845.598   6199.300   6024.433   5492.880   5475.516
## [31]   5071.906   4194.408   4169.282   3688.478   3611.481   3376.238
## [37]   3365.972   3203.616   3192.859   3179.853   3072.833   3014.386
## [43]   2996.338   2987.268   2984.135   2906.750   2672.857   2504.003
## [49]   2419.366   2386.851

plot(tsSSA$sigma, type="l")

Видим 13 “ненулевых” собственных числа.

rk = 13

par = parestimate(tsSSA, groups = list(1:rk), method = "esprit")

modulusReal = par$moduli
periodsReal = par$periods
o = order(abs(periodsReal), decreasing = TRUE)
r0 = reconstruct(tsSSA, groups = list(signal = 1:rk))$signal


len <- rk
vars <- matrix(nrow = len, ncol = rk)

for (i in 1:rk) {
  if (periodsReal[i] == Inf)
    vars[, i] <- modulusReal[i]^(1:len)
  else if (periodsReal[i] == 2)
    vars[, i] <- (-modulusReal[i])^(1:len)
  else if (periodsReal[i] > 0)
    vars[, i] <- 
      modulusReal[i]^(1:len) * sin(2 * pi * (1:len) / periodsReal[i])
  else
    vars[, i] <- 
      modulusReal[i]^(1:len) * cos(2 * pi * (1:len) / periodsReal[i])
}
lm0 <- lm(r0[1:len] ~ 0 + ., data = data.frame(vars))
coefs0 <- coef(lm0)
print(round(coefs0), digits = rk)

##        X1        X2        X3        X4        X5        X6        X7        X8 
##  -4965469  -5555248   7213618  54636156 -48963499      1651       739       205 
##        X9       X10       X11       X12       X13 
##       -65   1225627   -117262        NA        NA

idx <- seq(1, rk)
coefs.c.phase <- numeric(length(idx))
phases.c <- numeric(length(idx))
periods.c.phase <- numeric(length(idx))
moduli.c.phase <- numeric(length(idx))
for (i in seq_along(idx)) {
  
  periods.c.phase[i] <- periodsReal[idx[i]]
  moduli.c.phase[i] <- modulusReal[idx[i]]
  coefs.c.phase[i] <- coefs0[idx[i]]
  phases.c[i] <- atan2(coefs0[idx[i] + 1], coefs0[idx[i]])
  if(i == rk){
    phases.c[i] <- atan2(coefs0[idx[i]], coefs0[idx[i]])
  }
}

s = data.frame(periods = periods.c.phase, phases = phases.c, 
                 coefficients = coefs.c.phase, 
                 moduli =moduli.c.phase)
print(s)

##        periods      phases  coefficients    moduli
## 1    62.059089 -2.30019430 -4.965469e+06 1.0102164
## 2   -62.059089  2.22704139 -5.555248e+06 1.0102164
## 3   207.578092  1.43952547  7.213618e+06 1.0031497
## 4  -207.578092 -0.73069723  5.463616e+07 1.0031497
## 5          Inf  3.14155894 -4.896350e+07 1.0007054
## 6    11.995011  0.42101827  1.650952e+03 0.9986643
## 7   -11.995011  0.27033942  7.392869e+02 0.9986643
## 8     6.007629 -0.30780195  2.048739e+02 0.9982437
## 9    -6.007629  1.57084947 -6.513055e+01 0.9982437
## 10   38.009961 -0.09538464  1.225627e+06 0.9969873
## 11  -38.009961          NA -1.172619e+05 0.9969873
## 12   50.401313          NA            NA 0.9907013
## 13  -50.401313          NA            NA 0.9907013

ARIMA

Автоковариационная функция ACF

ts1 <- read.csv(file = "ts1.txt", header = TRUE, as.is = FALSE)
ts5 <- read.csv(file = "ts5.txt", header = TRUE, as.is = FALSE)

plot(ts1$x,  type = "l")

plot(ts5$x,  type = "l")

acf(ts1, main = "Функция для данных ts1")

acf(ts5, main = "Функция для данных ts5")

Функция PACF

По количеству столбцов, значимо отличающихся от нуля, можно определить количество параметров в авторегрессии.

pacf(ts1, main = "Функция для данных ts1")

pacf(ts5, main = "Функция для данных ts5")

Оценка параметров модели

fit1 <- auto.arima(ts1)
fit5 <- auto.arima(ts5)


fit1

## Series: ts1 
## ARIMA(2,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ar2     ma1     ma2     mean
##       -0.2652  0.6067  1.1344  0.1612  99.9179
## s.e.   0.0481  0.0335  0.0542  0.0454   0.1132
## 
## sigma^2 = 1.064:  log likelihood = -1448.07
## AIC=2908.15   AICc=2908.23   BIC=2937.59

fit5

## Series: ts5 
## ARIMA(3,0,0) with zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3
##       0.3103  -0.1605  0.0984
## s.e.  0.0222   0.0230  0.0223
## 
## sigma^2 = 0.3961:  log likelihood = -1910.39
## AIC=3828.78   AICc=3828.8   BIC=3851.18

forecastedValues1 <- forecast(fit1, 20)
forecastedValues5 <- forecast(fit5, 40)

plot(forecastedValues1, main = "Graph with forecasting FIT1",
col.main = "darkgreen")

plot(forecastedValues5, main = "Graph with forecasting FIT5",
col.main = "darkgreen")

ARIMA на реальных данных + (ARIMA VS SSA VS ETS)

fitArima <- auto.arima(unempl.male)
checkresiduals(fitArima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,2)(2,1,1)[12] with drift
## Q* = 20.36, df = 17, p-value = 0.2562
## 
## Model df: 7.   Total lags used: 24

summary(fitArima)

## Series: unempl.male 
## ARIMA(2,0,2)(2,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     sar1    sar2     sma1   drift
##       1.6149  -0.6340  -0.4319  0.1149  -0.0108  0.1072  -0.8538  5.8252
## s.e.  0.1395   0.1366   0.1483  0.0671   0.1022  0.0932   0.0909  3.6365
## 
## sigma^2 = 16160:  log likelihood = -2484.99
## AIC=4987.98   AICc=4988.44   BIC=5023.81
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.6993681 123.9684 91.59157 -0.34855 5.249927 0.1989636
##                      ACF1
## Training set 0.0005744349

fitArimaForecast <- forecast(fitArima, h = 24)
autoplot(fitArimaForecast)

tsSSABforecast <- bforecast(tsSSA, groups = list(c(1:4, 7:10)), len = 24, R = 50, only.new = F)


tsSSAETS <- hw(ts(unempl.male[1:(length(unempl.male)-24)], start=c(1948, 1), frequency=12), seasonal="multiplicative", h = 24)

autoplot(tsSSABforecast)

plot(tsSSAETS)

RMSE SSA vs ARIMA vs ETS

sd(unempl.male)

## [1] 794.9778

tsTest = tail(unempl.male, 24)

tsSSABforecast <- bforecast(tsSSA, groups = list(c(1:4, 7:10)), len = 24, R = 50, only.new = T)

mseSSA = mean((data.frame(tsSSABforecast)$Value - c(tsTest))**2)

mseARIMA = mean((c(fitArimaForecast$mean) - c(tsTest))**2)

mseETS = mean((c(tsSSAETS$mean) - c(tsTest))**2)

RMSE SSA

print(sqrt(mseSSA))

## [1] 1399.446

RMSE ARIMA

print(sqrt(mseARIMA))

## [1] 1024.97

RMSE ETS

print(sqrt(mseETS))

## [1] 889.6957

SSA и разладка

library(Rssa)
library("lattice")
firstHalf = c(10*sin(2*pi*(1:72)/6))
secondHalf = c(30*sin(2*pi*(73:144)/6))
tsdestruct = c(firstHalf,secondHalf)
plot(tsdestruct, type = "l", col = "red")

s <- ssa(tsdestruct, L = 12)
w <- wcor(s, groups = 1:10)
plot(w)

plot(s, type = "vectors", idx=1:10)

r <- reconstruct(s, groups = list(Trend = c(1)))

tsdestruct_res <- residuals(r)
N <- length(tsdestruct_res)
rank <- 2
periods <- function(M, L) {
  ts(sapply(1:(N - M),
            function (i) {
              s <- ssa(tsdestruct_res[i:(i + M - 1)], L = L)
              par <- parestimate(s, groups = list(c(1:rank)), 
                                 method = "esprit")
              abs(par$periods[1])
            }),
     start = time(tsdestruct)[M + 1], delta = 1)
}
per12 <- periods(12, 6)
per24 <- periods(24, 12)
lattice::xyplot(plot.method = "xyplot",per12, type = "l")

M <- 12; L <- M / 2
hm <- hmatr(tsdestruct_res, B = M, T = M, L = L, neig = rank)
plot(hm)

Timeseries

DM

2023-11-07

Библиотеки

Периодограммы

Эффект растекания синуса

Генерация белого шума

Генерация красного шума

Построение периодограммы для реальных данных

График безработицы в США

Периодограмма безработицы в США

Фильтры

Функция для фильтра скользящего среднего

АЧХ

Применение фильтра скользящего среднего к модельным данным

Задаем модель

Исходная периодограмма

Применяем фильтр

Периодограмма MA(5)

Периодограмма MA(10)

АЧХ фильтра скользящего среднего

Применение фильтра к реальным данным

SSA

Выделение тренда с помощью SSA

Матрица взвешенных корреляций

Собственные вектора

FOSSA

FOSSA: матрица корреляций

FOSSA: собственные вектора

FOSSA: итоговое разложение

Корреляции SSA vs. FOSSA

Сезонности SSA vs. FOSSA

Тренд SSA vs. FOSSA

Различные способы выделения тренда

Тренд по POLY

Тренд по LOESS

Тренд по HP

Итоговый график со всеми способами выделения тренда

Сравнение методов разложения

STL

Classical seasonal decomposition

Итоговый график трендов разных разложений

Примеры применения Toeplitz SSA и SSA with projections

SSA with projections

Toeplitz SSA

Прогнозы: рекурентный и векторный

Рекурентный прогноз - rforecast

C отсечением периода

Векторный прогноз - vforecast

C отсечением периода

Восстановление ряда по ЛРФ

Ряд 1: Экспонента

Ряд 2: Линейная функция

Восстановление ЛРФ реального временного ряда

ARIMA

Автоковариационная функция ACF

Функция PACF

Оценка параметров модели

ARIMA на реальных данных + (ARIMA VS SSA VS ETS)

RMSE SSA vs ARIMA vs ETS

RMSE SSA

RMSE ARIMA

RMSE ETS

SSA и разладка